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i Abstract. We study numerically the late-time behaviour of the coupled Einstein 

Yang-Mills system. We restrict ourselves to spherical symmetry and employ Bondi- 
' like coordinates with radial compactification. Numerical results exhibit tails with 

exponents close to —4 at timelike infinity and —2 at future null infinity . 
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1. Introduction 



(N 

OO , Radiating systems relax to equilibrium by dissipating energy to infinity. The fall-off 
^ ! properties of the field at late times are governed by so-called radiation tails. These tails 
emerge from primary outgoing radiation that is backscattered. This far-field effect is 
either due to a background or an effective potential produced by the nonlinearities of 
OO ! the radiation field itself, or in general, a mixture of both. 

The classical fall-off properties were based on liner perturbations about a given 
background [1]. However, Bizon [2] has pointed out, based on previous results [31 HI [5], 
^ i that for certain nonlinear systems the nonlinear tails may dominate the long-time 
5^ I behaviour, i.e. these tails fall off more slowly in time than the linear perturbations. As an 
example, he studied the spherically symmetric Yang-Mills equations on Minkowski and 
Schwarzschild spacetimes [6]. While linear perturbation theory predicts a power law 
decay due the backscattering off the Schwarzschild curvature, the nonlinear part decays 
only as for observers near timelike infinity. This slower fall-off was also observed 
numerically. More recently, these results were confirmed and extended to the late-time 
behaviour at future null infinity by Zenginoglu [7], showing that tails die off as on 
Schwarzschild spacetime. 

These calculations are on a given background and the question arises what happens 
for the coupled Einstein- Yang-Mills system. In the fully coupled case it is difficult to 
disentangle the different contributions. If one writes down a perturbation expansion 
starting from flat spacetime, then the flrst order perturbation of the YM fleld evolves 
in an effective potential produced by the back reaction of the YM-fleld to the metric. 
In addition, there will be the nonlinear effects from the YM equation already present in 



Tails for the Einstein- Yang-Mills system 



2 



flat space. In this sense all tails are nonlinear, however, it is not evident what the overall 
fall-off behaviour is. The aim of this paper is to answer this question numerically. 

It is known that the possible stable endstates of (spherically symmetric) collapse 
for the Einstein- YM system are either the formation of a black hole or dispersion to 
flat space. Black holes come in two types: pure black holes, i.e. with all of the YM 
field radiated away or coloured pj. In addition, there exist soliton solutions found by 
Barnik and McKinnon However, both solitons and colored black holes turned out 
to be unstable under linear perturbations. For technical reasons we restrict ourselves 
to subcritical evolutions which do not form black holes. However, from the results in 
P, [^, we expect our findings to hold also for generic spherically symmetric initial data 
when a black hole is formed. 

We tackle this problem by numerically solving a characteristic initial value problem 
with constraints and employ radial compactification to allow investigation of tails at 
future null infinity. First we present the model, derive the equations of motion and 
write the Yang- Mills wave equation as an advection equation plus a constraint. Second, 
we discuss the numerical methods used for solving the given system and check the 
convergence of the code. Finally, we present results for late-time tails on Minkowski 
background and for the coupled Einstein- Yang- Mills system. 



2. Model 



We assume spherical symmetry with a regular center and choose Bondi-like coordinates 
{u, r, 9, 0} based upon outgoing null hypersurfaces u = const with the line-element 



r 



(1) 



We consider the Yang-Mills theory with the gauge group SU(2) and assume the 
magnetic ansatz for the gauge connection [6l [Tl] 

A = wT^de + (cot Ot' + wr"^) sin ed(j), (2) 

where w = w{u,r) is the Yang-Mills field and the t"" are the spherical generators of 
su{2), normalized such that [t'^,t^] = ie"'"^r'^, where a,b,c ^ {r,6,(j)}. They are related 
to the Pauli matrices a" via r" = cr"/2. 

The Yang-Mills field strength (or curvature), F = dA + A /\ A, then becomes 



F = {wdu + w'dr) A [r^dO + r"^ sin 



(1 - w^) r"rf0 Asin^rf0, (3) 



where w and w' denote partial derivatives of w{u, r) with respect to u and r, respectively. 
Using the above ansatz the trace of the Yang-Mills curvature becomes 



tr (F^^F, 



-2/3 



2ww' - ^iw'f 
r 



+ 



w 



2\2 



(4) 



where Greek indices range over the four spacetime dimensions and Latin indices are 
group indices. The action for the Yang-Mills field coupled to Einstein's equations is [TT] 



S 



d Xy 



R 
IQttG 



l2 



(5) 
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The Yang-Mills wave equation is obtained by varying the action with respect to the 
Yang-Mills field w 



V 



V 



-2w' + [ — ] w' + —w" + -^w{l 



w 



0, 



(6) 



while variation with respect to the metric functions, f3 and V, yields two constraint 
equations 



v 



w 



2\2 



SvrG iw 



i\2 



(7) 

r 

The coupling constant [G/e^] has dimension of length'^. Since it is not 
dimensionless, changing the coupling constant does not give rise to a one-parameter 
family of theories, but only changes the scale. To simplify the equations, we choose 

5? ^ 1. (9) 



The final form of the hypersurface equations then becomes 



V 



1 



l\2 



(10) 

(11) 



The regularity condition to be imposed on the Yang-Mills field at the origin is 

w = ±l + 0{r'^) (12) 

while the gauge {u is chosen to be proper time at the regular center) and regularity 
conditions on the metric functions are 

/? = 0(r2) (13) 
V/r = l + 0{r'^). (14) 

For the study of tail behaviour it is crucial to introduce a new field variable 

w := w — 1 (15) 

as a perturbation of one of the Yang-Mills vacua w = ±1, (i.e. where the field 
strength vanishes), in order to avoid the tails being swamped by accumulated numerical 
errors. Since the matter field equation and the constraint equations are invariant under 
refiection symmetry, w{u, r) — > w{u, r) := —w{u, r), we may specialize to one of the two 
vacua without loss of generality. 

There are a number of different ways to go about solving the Yang-Mills wave 
equation Qj in Bondi coordinates, e.g. the diamond integral approach due to Gomez 
and Winicour [121 [13] which we used in jlO] or Goldwirth and Piran [13] and Garfinkle's 
[T3] way of rewriting the equation with a total time derivative along the ingoing null 
geodesies. The latter allows for employing standard method of lines (MOL) techniques, 
i.e. one first discretizes the spatial derivatives (albeit non-equispaced) , which in turn. 
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yields a system of coupled ordinary differential equations (ODEs) that can be solved 
using a standard ODE solver. In this approach to solving the characteristic initial value 
problem, the gridpoints usually move along the ingoing null geodesies, which implies a 
nontrivial origin treatment, with Taylor expansions for increased accuracy. 

In contrast to the methods mentioned above, we prefer to simply introduce a new 
evolution variable 

h := w^r (16) 

so as to eliminate the mixed ur derivative in equation ([6]). In addition, we also keep the 
locations of gridpoints (in time) at fixed values of r. Here, MOL discretizations, using 
standard stencils for equidistant grids, are applicable and the treatment of the origin is 
trivial modulo boundary conditions. Moreover, we are here not interested in strong field 
regions, where the focussing of ingoing null geodesies would provide a natural increase 
of resolution near the center. Rather, we want to study late time tails for subcritical 
evolutions, which entails tracking the field at locations of constant r through time. 
The evolution equation then becomes 



where 



F{h) = 2h + 3h^ + h^, 



and h = w. 

We now have an added constraint to solve (in addition to the two geometry 
equations (fTOD. ffTT]) ) 

h = [ h{u,r)dr. (19) 
Jo 

Note that equation (|T71) is of advection type. For flat space and without the YM 
self-interaction term, it reduces to 

h = h', (20) 
which is equivalent to the fiat space wave equation for 

1 r 

= - / h(u,r)dr. (21) 
r Jo 

Given initial data h{u = 0, r) = /(r) with r G M the solution of the advection equation 
(l20l) is simply 

h{u,r) = f{r + ^u), reR,u>0. (22) 

The characteristic curves r + u/2 = const are purely ingoing and thus, there is no 
boundary condition at the origin of spherical symmetry. The outgoing characteristic of 
the wave equation comes into play via the constraint equation ( 12T|) . 

The outer boundary needs a different treatment. While it is possible to specify 
an outgoing wave boundary condition we prefer to use compactification, which has the 
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advantage of being able to observe late-time behaviour at future null infinity. As in [10] 
we introduce a compactified radial coordinate 

r 



X 



1+r' 



(23) 



which maps r G [0, oo] i-^ a; G [0, 1]. In addition, we use the Misner-Sharp mass-function 



1 - ^e-^-^ 

r 



(24) 



as an evolution variable, thereby eliminating V. This is necessary, since the compactified 
constraint equation for V is singular at future null infinity J^~^, similar to the scalar 
field case treated in llOi. Instead of h we introduce a new field variable 



h := w^^ = h^^, 

which finally leads to a manifestly non-singular evolution system. 
The evolution equation then becomes 



(25) 



,2/3 



2m- 



X 



X 



'l~xYh 



.2/3 



x^ 



and the constraints are 

/ hdx 
Jo 

(1 -x)3 



h 



m . 



X 

1 - 2m 



1 — X 

X 



{i-x)\hY + 



2a;2 



(26) 

(27) 
(28) 
(29) 



The regularity conditions at the origin for the compactified scheme become, using 
r = 0(x), 

h = 0{x) h = 0{x^) (30) 
(3 = O(x^) m = 0{x^) (31) 

There is no boundary condition at future null infinity, since it is a characteristic. 



3. Numerics 



From a numerical point of view, we do not even need to know that we are solving a 
characteristic initial value problem. Rather, we may simply take the evolution system 
and solve the nonlinear advection equation (126|) using a standard MOL approach in 
conjunction with computing the constraint ODEs ( |271) . (!28i) . and (l29l) . 

We may discretize the advection term in fl26l) by centered, fully upwind or upwind- 
biased stencils. Combined with an ODE integrator for time, such as, the classical 4th 
order Runge-Kutta method (RK4), the resulting schemes will then exhibit different 
numerical errors. One has to make a choice between added dispersion, in the case of 
centered approximations, and added dissipation for upwind schemes. These errors affect 



Tails for the Einstein- Yang-Mills system 



6 



mostly high frequency components of the solution and can be minimized by using high 
order approximations. 

In the interior of the grid we have chosen a 6th order upwind-biased scheme with 
the stencil 

dH 2ifj_2 — 24i7i_i — 35i/i + 80i/i+i — 30ifj+2 + 8i/j+3 — , a-. 

a^^"^) = GOAx + ^'^^ 

where Hi = H{xi). The above upwind stencil is the one closest to the centered stencil 
and it also has the lowest error term among the upwind stencils [12] . The two alternate 
upwind stencils do not lead to stable evolutions in our case. The weights for such 
finite-difference stencils can conveniently be computed in a computer algebra package, 
such as Mathematica, using Fornberg's compact algorithm [17]. In contrast to centered 
schemes, where it is often necessary to add some artificial dissipation to have numerical 
stability, the dissipation is already "built-in" in our chosen scheme. 

The term F(h)/x^ in fl26|l forces a regularity boundary condition at the origin, so 

that 

/i(u,a; = 0) = 0. (33) 

We enforce it by choosing initial data that satisfy this condition (to machine precision) 
and then make sure that its time derivative is zero, i.e. hu{u, x = 0) = for all timesteps. 
We choose Gaussian initial data 

/i(0, x) = A exp [-200(2; - 1/2)2] _ (^34^) 

The constraint equations for h and (3 are integrated using cumulative Newton- 
Cotes quadrature rules of order 6, which are given in the appendix. Since the right 
hand side of the constraint equation for the Misner-Sharp mass-function (1291) depends 
on the unknown m, we integrate it with RK4 and use 4th order polynomial interpolation 
for computing h and h at points Xj+1/2 in between actual gridpoints, as required by the 
Runge-Kutta method. 

In comparison to the established methods for solving characteristic initial value 
problems [T2l [T3l HH |15|, the method used here has a number of advantages for the 
problem at hand. It relies on a standard MOL discretization using equidistant stencils 
and thus allows for the use of high order schemes. Moreover, the boundary treatment 
is simple and does not require Taylor series near the origin. It also allows us to track 
the field at lines of constant r without further interpolation. 

In general, we expect the code to be 4th order convergent, see figure [TJ For small 
data, however, max2m/r is also small. Since m only appears in the wave equation in 
this combination, errors in the computation of m may be allowed to be larger than those 
in other fields, and we may therefore have close to 6th order spatial accuracy, in this 
case. If the Courant number 

(-) 

is small enough, so that the errors from the RK4 integrator are of order 0(Aa;^), then, 
we may in fact achieve 6th order convergence. It is, however, impractical to have small 
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Courant numbers for very long time evolutions. Clearly, evolutions using high Courant 
numbers complete faster and need less timesteps than evolutions using lower Courant 
numbers. Therefore, accumulation errors should also be somewhat less severe for high 
Courant numbers. Moreover, the classical RK4 solver exhibits additional damping [18] of 
high frequency modes near the stability limit, which is beneficial to numerical stability. 
For weak fields, the Yang-Mills advection equation fl26l) approximately reduces to 

k = ^{l-xfh^. (36) 

After freezing the nonconstant coefficient (1 — x)^ at its maximum 1, Fourier analysis 
leads to the stability limit (similar to the 6th order centered case discussed in |19] ) 



C=-^< 1.29 (37) 

for a MOL discretization with the 6th order upwind-biased stencil (1321) and RK4 as 
time-integrator. For the coupled system, we have found that Courant numbers of up to 
C < 1.25 yield stable evolutions. 



7 




Figure 1. This figure shows the convergence factor C/ = log2 |j,-2ooo feiooo|| 
Z2-norni of the solution for the coupled EYM-system with grid resolutions of 500, 1000 
and 2000 gridpoints and a Courant number of C = 1.25. Initially the code is 4th 
order convergent. Since the initial amplitude A ~ 0.28 is quite large and 2m/ r w 0.6 
the spatial accuracy is only 4th order. Later in the evolution, as the field disperses, 
2m /r becomes very small and the spatial discretization becomes roughly 6th order 
accurate. At even later times the convergence order decreases slowly - probably due 
to the standard degrading of convergence in long evolutions. 

We have used numerical Python to conveniently automate the determination of tail 
exponents via fitting, while the core of the code was written in C-I--I-. 
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4. Results 

Perturbation theory predicts the late-time behaviour of the solution to be 

h{u,x)^CuP, (38) 

where p is the tail exponent. A first test case for our code is to correctly reproduce 
the know tail behaviour on Minkowski background [6] . In figure [2] we find that the 
tail exponent tends to p = —4 for observers near timelike infinity and the exponent 
p = —2 for future null infinity This also corresponds to the decay found on 

Schwarzschild backgrounds [7J. In terms of the compactified radial coordinate x, the 
observers are located at x = (1, 0.9995, 0.999, 0.995, 0.99, 0.95, 0.9, 0.8, 0.7, 0.5). For the 
results presented here, we have used 10000 spatial gridpoints, a Courant number of 
C = 1.25 and an initial data amplitude of A = 0.01. For bigger, but still subcritical 
(in the coupled case), amplitudes and/or smaller Courant numbers, the tail decay is 
essentially the same. 

Figure [3] encodes our main results showing the late-time behavior for the coupled 
Einstein- Yang-Mills system. We find essentially the same fall-off as for the Yang- 
Mills field on Minkowski background. As mentioned in the introduction, in terms of 
a perturbative approach, tails are generated on the one hand by the the nonlinearity 
of the YM field itself, and, on the other by the contribution of the field to the metric. 
What we see numerically is a superposition of both effects which we can not separate. 
Hod [20] has studied linear wave tails in time dependent potentials. He finds that for a 
certain class of potentials that go to zero asymptotically in time, the fall-off behaviour 
of the tails is a power law depending on the time dependence of the potential. 

5. Conclusion 

Using Bondi-like coordinates and radial compactification we have written the Einstein- 
Yang-Mills system as an advection equation plus three constraints. In this form, MOL 
discretizations are straightforward to apply, the origin treatment is easy and the outer 
boundary, J^'^, being a characteristic does not require boundary data. Compared to 
the diamond integral scheme [121 [13] Goldwirth, Piran and Garfinkle's [HI [15] way 
of solving characteristic initial value problems, the approach used here is very clean, 
simple to implement and allows the use of high order schemes. 

We have found that the spherically symmetric coupled Einstein- Yang-Mills system 
shows the same fall-off behaviour at late times as Yang-Mills on Minkowski or 
Schwarzschild backgrounds. Although such a result could have been expected it is 
by no means evident, because so far it is not known how tails arising from the back 
reaction of the the Yang-Mills field to the metric decay. Our results indicate that they 
decay as fast or faster than the nonlinear tails on Minkowski background. 

Our compactified code has allowed us to also study the fall-off behavior at future 
null infinity. As for the coupled Einstein massless scalar field [10], the fall-off on ^"""is 
slower than for observers approaching timelike infinity. Since realistic observers are 
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Figure 2. The upper plot shows the decay of the field h at x — const versus retarded 
time u, while the lower plot depicts the respective tail exponents for the same evolution 
on Minkowski background. On ^+the tail exponent is close to p = —2. Observers 
located at finite x (or r) approach timelike infinity i+ with the exponent p = —4 for 
late times. Observers closer to the center approach this value faster than those closer 
to J^+. 



located only at finite distances from the center, what then is the practical relevance to 
know the decay conditions on ^'*"? It has been pointed out in [10] and also in [7], that 
for astrophysical observers, the relevant decay rate is the one along null infinity. This 
has to do with the observation that the tail exponents for observers far out start close 
to the exponent on ^^ajidi. only slowly decrease to the value for timelike observers. 
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Figure 3. Similar to figure [2] tliese plots depict tlie late-time tails and the respective 
exponents for the coupled Einstein- Yang-Mills system. The exponents coincide with 
the behaviour on Minkowski space, being p — —2 at .y^andp — —4 at , respectively. 
The zero-crossing in /i at u « 50 depends on the initial data amplitude, i.e. the field 
goes through zero earlier for smaller amplitudes. 
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Appendix 

The quadrature formulas below have been obtained by simply integrating the (quartic) 

interpolating polynomial P(/|a;j_i, a;,, Xj+i, a;i+2, Xj+s) on an equidistant grid with 
spacing h over the intervals [xi-i,Xi] until [a;i+2, Xj+s], respectively. 

£ fdx = A (251/,_i + 646/, - 264/,+! + 106/,+2 - 19/^+3) + 0{h') (A.l) 

£ ^ /c^^ = ^ (-19/^-1 + 346/, + 456/,+i - 74/,+2 + ll/^+a) + 0{h') (A.2) 

fdx = A _ 74/^ + 456/.+, + 346/,+2 - m+s) + 0(/i^) (A.3) 



/^^^ = ^ (-19/i-i + 106/, - 264/,+i + 646/,+2 + 251/,+3) + 0(/i^) (A.4) 

Summing these formulas together, yields the classical Boole's or Milne's rule. 
/ fdx = — (7/,_i + 32/, + 12/,+i + 32/,+2 + 7/,+3) + 0(/i^) (A.5) 
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